### Load and Merge Data ###
# rm(list = ls())
# setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
warnings("off")
library(dplyr)
library(lfe)

df_main <- fread("output/ConstraintMeasures_A.csv")
df_TFP <- fread("main_cleaned.csv")
df_marginQ <- haven::read_dta("data/marginalq_manufacturing_73_19.dta")
df_IntanQ <- fread("data/IntanQ_Peters_Taylors.csv")

# winsorize intangible Q. Note that the KZ and WW indexes have the same winsorization step in the pre-processing step in file "FinancialConstraintMeasures_Annual.R". 
df_IntanQ <- df_IntanQ %>% 
  mutate(q_tot = wins(q_tot))

df_main <- df_main %>% 
  mutate(gvkey = as.numeric(gvkey)) %>% 
  left_join(df_TFP %>% select(gvkey, fyear, TFP), by = c("gvkey", "year"="fyear"))%>% 
  left_join(df_marginQ, by = c("gvkey", "year"="fyear"))%>% 
  left_join(df_IntanQ %>% select(gvkey, fyear, q_tot), by = c("gvkey", "year"="fyear"))

df_main <- df_main %>% 
  filter(!is.na(WW),
         !is.na(KZ),
         !is.na(TFP))

# residuzlize 1 (average Q)
lm_WWRes1 <- felm(WW~TobinQ|gvkey+year|0|gvkey+year, df_main)
lm_KZRes1 <- felm(KZ~TobinQ|gvkey+year|0|gvkey+year, df_main)
df_main <- df_main %>% mutate(WW_res1 = lm_WWRes1$residuals)
df_main <- df_main %>% mutate(KZ_res1 = lm_KZRes1$residuals)

# residuzlize 2 (marginal Q)
# Note: Marginal Q has a lot of missing values
tmp <- df_main %>% 
  filter(!is.na(marginalQ))

lm_WWRes2 <- felm(WW~marginalQ|gvkey+year|0|gvkey+year, tmp)
lm_KZRes2 <- felm(KZ~marginalQ|gvkey+year|0|gvkey+year, tmp)
tmp <- tmp %>% mutate(WW_res2 = lm_WWRes2$residuals)
tmp <- tmp %>% mutate(KZ_res2 = lm_KZRes2$residuals)

df_main <- df_main %>% 
  left_join(tmp %>% select(gvkey, year, WW_res2, KZ_res2),
            by = c("gvkey", "year"))

# residuzlize 3: use total Q from Peters and Taylor (2016)
tmp <- df_main %>% 
  filter(!is.na(q_tot))

lm_WWRes3 <- felm(WW~q_tot|gvkey+year|0|gvkey+year, tmp)
lm_KZRes3 <- felm(KZ~q_tot|gvkey+year|0|gvkey+year, tmp)
tmp <- tmp %>% mutate(WW_res3 = lm_WWRes3$residuals)
tmp <- tmp %>% mutate(KZ_res3 = lm_KZRes3$residuals)

df_main <- df_main %>% 
  left_join(tmp %>% select(gvkey, year, WW_res3, KZ_res3),
            by = c("gvkey", "year"))


df_main <- df_main %>% 
  mutate(financial_capacity_KZ_without_Q = -KZ_withoutQ)
df_main <- df_main %>% 
  mutate(financial_capacity_KZ_res1 = -KZ_res1 )
df_main <- df_main %>% 
  mutate(financial_capacity_KZ_res2 = -KZ_res2 )
df_main <- df_main %>% 
  mutate(financial_capacity_KZ_res3 = -KZ_res3 )
df_main <- df_main %>% 
  mutate(financial_capacity_WW_res1 = -WW_res1 )
df_main <- df_main %>% 
  mutate(financial_capacity_WW_res2 = -WW_res2 )
df_main <- df_main %>% 
  mutate(financial_capacity_WW_res3 = -WW_res3 )


### Regression Table ###
lm1 <- felm(financial_capacity_KZ_res1~TFP|gvkey+year|0|gvkey+year, df_main)
lm2 <- felm(financial_capacity_KZ_res2~TFP|gvkey+year|0|gvkey+year, df_main)
lm3 <- felm(financial_capacity_KZ_res3~TFP|gvkey+year|0|gvkey+year, df_main)
lm4 <- felm(financial_capacity_WW_res1~TFP|gvkey+year|0|gvkey+year, df_main)
lm5 <- felm(financial_capacity_WW_res2~TFP|gvkey+year|0|gvkey+year, df_main)
lm6 <- felm(financial_capacity_WW_res3~TFP|gvkey+year|0|gvkey+year, df_main)
# lm5 <- felm(financial_capacity_KZ_without_Q~TFP|gvkey+year|0|gvkey+year, df_main)
lm.l <- list(lm1,lm2,lm3,lm4,lm5,lm6)
stargazer(lm.l, type = "text", omit.stat = c("F", "ser", "adj.rsq"), dep.var.caption = "Dep var: Financial capacity index", dep.var.labels = rep(c("average Q", "marginal Q", "total Q"), 2) )
stargazer(lm.l, omit.stat = c("F", "ser", "adj.rsq"), dep.var.caption = "Dep var: Financial capacity index", dep.var.labels = rep(c("average Q", "marginal Q", "total Q"), 2), add.lines = list(c("Firm Fixed Effect", rep("Yes", 6)), c("Time Fixed Effect", rep("Yes", 6))),   out = "../Tables/TableA1_regression_results.txt")



### Plot of financial capacity index ###
figure_width = 4
figure_height = 3
warnings("off")
df_main %>% 
  mutate(Y = financial_capacity_KZ_res1) %>% 
  filter(!is.na(Y),
         !is.na(TFP)) %>% 
  mutate(bin = as.character(ntile(TFP, n=10))) %>%
  group_by(bin) %>%
  summarise(n = n(),
            FC_mean = mean(Y), 
            FC_sd = sd(Y), 
            TFP_mean = mean(TFP)) %>%
  mutate(FC_se = FC_sd/sqrt(n)) %>%
  ggplot(aes(x=TFP_mean, y=FC_mean)) + 
  geom_point()+
  geom_errorbar(aes(ymax=FC_mean+2*FC_se, ymin=FC_mean-2*FC_se))+
  scale_x_continuous(name = "TFP") +
  scale_y_continuous(name = "financial capacity") +
  theme_bw()
ggsave( paste0(FigurePath, "/FigureA6a_TFP_KZres1.pdf"), width=figure_width, height=figure_height)

df_main %>% 
  mutate(Y = financial_capacity_KZ_res2) %>% 
  filter(!is.na(Y),
         !is.na(TFP)) %>% 
  mutate(bin = as.character(ntile(TFP, n=10))) %>%
  group_by(bin) %>%
  summarise(n = n(),
            FC_mean = mean(Y), 
            FC_sd = sd(Y), 
            TFP_mean = mean(TFP)) %>%
  mutate(FC_se = FC_sd/sqrt(n)) %>%
  ggplot(aes(x=TFP_mean, y=FC_mean)) + 
  geom_point()+
  geom_errorbar(aes(ymax=FC_mean+2*FC_se, ymin=FC_mean-2*FC_se))+
  scale_x_continuous(name = "TFP") +
  scale_y_continuous(name = "financial capacity") +
  theme_bw()
ggsave( paste0(FigurePath, "/FigureA6e_TFP_KZres2.pdf"), width=figure_width, height=figure_height)


df_main %>% 
  mutate(Y = financial_capacity_KZ_res3) %>% 
  filter(!is.na(Y),
         !is.na(TFP)) %>% 
  mutate(bin = as.character(ntile(TFP, n=10))) %>%
  group_by(bin) %>%
  summarise(n = n(),
            FC_mean = mean(Y), 
            FC_sd = sd(Y), 
            TFP_mean = mean(TFP)) %>%
  mutate(FC_se = FC_sd/sqrt(n)) %>%
  ggplot(aes(x=TFP_mean, y=FC_mean)) + 
  geom_point()+
  geom_errorbar(aes(ymax=FC_mean+2*FC_se, ymin=FC_mean-2*FC_se))+
  scale_x_continuous(name = "TFP") +
  scale_y_continuous(name = "financial capacity") +
  theme_bw()
ggsave(paste0(FigurePath, "/FigureA6c_TFP_KZres3.pdf"), width=figure_width, height=figure_height)

df_main %>% 
  mutate(Y = financial_capacity_WW_res1) %>% 
  filter(!is.na(Y),
         !is.na(TFP)) %>% 
  mutate(bin = as.character(ntile(TFP, n=10))) %>%
  group_by(bin) %>%
  summarise(n = n(),
            FC_mean = mean(Y), 
            FC_sd = sd(Y), 
            TFP_mean = mean(TFP)) %>%
  mutate(FC_se = FC_sd/sqrt(n)) %>%
  ggplot(aes(x=TFP_mean, y=FC_mean)) + 
  geom_point()+
  geom_errorbar(aes(ymax=FC_mean+2*FC_se, ymin=FC_mean-2*FC_se))+
  scale_x_continuous(name = "TFP") +
  scale_y_continuous(name = "financial capacity") +
  theme_bw()
ggsave(paste0(FigurePath, "/FigureA6b_TFP_WWres1.pdf"), width=figure_width, height=figure_height)

df_main %>% 
  mutate(Y = financial_capacity_WW_res2) %>% 
  filter(!is.na(Y),
         !is.na(TFP)) %>% 
  mutate(bin = as.character(ntile(TFP, n=10))) %>%
  group_by(bin) %>%
  summarise(n = n(),
            FC_mean = mean(Y), 
            FC_sd = sd(Y), 
            TFP_mean = mean(TFP)) %>%
  mutate(FC_se = FC_sd/sqrt(n)) %>%
  ggplot(aes(x=TFP_mean, y=FC_mean)) + 
  geom_point()+
  geom_errorbar(aes(ymax=FC_mean+2*FC_se, ymin=FC_mean-2*FC_se))+
  scale_x_continuous(name = "TFP") +
  scale_y_continuous(name = "financial capacity") +
  theme_bw()
ggsave(paste0(FigurePath, "/FigureA6f_TFP_WWres2.pdf"), width=figure_width, height=figure_height)

df_main %>% 
  mutate(Y = financial_capacity_WW_res3) %>% 
  filter(!is.na(Y),
         !is.na(TFP)) %>% 
  mutate(bin = as.character(ntile(TFP, n=10))) %>%
  group_by(bin) %>%
  summarise(n = n(),
            FC_mean = mean(Y), 
            FC_sd = sd(Y), 
            TFP_mean = mean(TFP)) %>%
  mutate(FC_se = FC_sd/sqrt(n)) %>%
  ggplot(aes(x=TFP_mean, y=FC_mean)) + 
  geom_point()+
  geom_errorbar(aes(ymax=FC_mean+2*FC_se, ymin=FC_mean-2*FC_se))+
  scale_x_continuous(name = "TFP") +
  scale_y_continuous(name = "financial capacity") +
  theme_bw()
ggsave(paste0(FigurePath, "/FigureA6d_TFP_WWres3.pdf"), width=figure_width, height=figure_height)


